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Abstract 

The algorithm for Dissipative Particle Dynamics (DPD), as modified by Espafiol and 
Warren, is used as a starting point for proving an if-theorem for the free energy and de- 
riving hydrodynamic equations. Equilibrium and transport properties of the DPD fiuid are 
exphcitly calculated in terms of the system parameters for the continuous time version of 
the model. 

PACS. 05.20 Dd Kinetic theory; 05.70 Ln Nonequilibrium thermodynamics, irreversible 
processes; 47.11+j Computational methods in fluid dynamics. 







In recent years several new simulation methods have been proposed for studying dynami- 
cal and rheological properties of complex fluids, on time scales which are difficult to reach by 
conventional molecular dynamics methods. These new techniques include lattice gas cellular 
automata (LGCA), the lattice Boltzmann equation (LBE), and dissipative particle dynamics 
(DPD). 

The goal of this letter is a theoretical analysis of the equilibrium and transport properties 
of a DPD system. This method was introduced by Hoogerbrugge and Koelman and 
modified by Espanol and Warren to ensure a proper thermal equilibrium state. There 
exist many recent applications of this new technique to concentrated colloidal suspensions, 
dilute polymer solutions and phase separation in binary mixtures ||^. Of course, to study the 
rheology of colloidal particles and polymers, suspended in a DPD fluid, one needs to model 
in addition the hydrodynamic interactions between these objects, transmitted through the 
fluid which is not the goal of this article. 

We briefly recall the basic elements of the simulation method. The DPD particle-based 
algorithm considers N point particles which model a fluid out of equilibrium and conserve 
mass and momentum, but not energy. Positions and velocities are continuous, but time is 
discrete and incremented in time steps 5t, as in LGCA and LBE. The algorithm consists of 
a collision step and a propagation step. In the collision step, the velocity of each particle 
is updated according to its interaction with particles inside a sphere of radius Rq through 
conservative, frictional and random forces. In the subsequent propagation step of fixed length 
6t all particles move freely to their new positions. 

The equations of motion for DPD with a finite step 6t are given by [0 

Vi{t + St) - v,{t) = ai(t, 6t) = J2 {^ij^t + t^ij^Wij] 

ri{t + 5t)-ri{t) = Vi{t + 6t)6t, (1) 

where Xi = {vi, Vi} is the phase of the i-th particle {i = 1,2, . . . , N). The interparticle force 
contains a systematic part k and a random part (t6W. The systematic force kij = Fij/m — 
has a conservative part ~ F and a dissipative part 7. The random force is described 
by a stochastic variable 6W, defined through 6Wij = ^1^^^ dT^ij^r), where iijit) is Gaussian 
white noise with an average = and a correlation ^ij(t)^ke(t') = {SikSje + 6ie6jk)S(t — t'), 
and 5{t) is a Dirac delta function. The forces are defined as Fij = —d(j){Rij)/dri,j^j = 
{vi — Vj) and (Tij = awR{Rij)Rij, where 0(-R) is a pair potential, wd{R) 
and wr{R) are positive weight functions, and 7 and a are respectively the strength of the 
friction and the noise. The weight function vanishes outside a finite range Rq (here chosen 
to be equal for wd and wr). Moreover, Rij = Vi — Vj, and a hat denotes a unit vector. 
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The interparticle forces kij, Fij, and cTij are antisymmetric under interchange of i and 
j, which guarantees the conservation of total momentum, which is a prerequisite for the 
existence of a slowly chainging local fluid velocity and the validity of the Navier Stokes 
equation. 

Formally, the DPD algoirthm defines a microscopic N particle system. However, the 
introduction of noise and dissipation represents a coarse grained description. Consequently, 
the "DPD particles" should not be interpreted as molecules, but as some mesoscopic degree 
of freedom of the fluid, referred to as "lumps." If to oc l/'^nR^ denotes the characteristic 
molecular time scale in DPD with n = N/V the number density, d the number of dimensions 
and 7 the friction constant, then to is considered to be large compared to any molecular time 
scale. In this description the dominant interactions are friction (~ 7) and random noise 
(~ cr), whereas the conservative forces can be interpreted as weak forces of relatively long 
range, which can be treated perturbatively. At sufficiently small temperature, they are 
responsible for phase transitions from liquid to crystalline order. Here they will be set equal 
to zero (7 large). The random forces act effectively as repulsive forces to prevent collapse of 
DPD particles. 

Define the single particle and pair distribution functions as 

f{x,t) = {J26ix-x,m; f^'\x,x',t) = {Y^6ix-x,it))6ix'-x,{t))), (2) 

where < ■ ■ ■ > represents an average over an initial ensemble of the A^-particle system. 
Combination of (|^) and then yields: 

f{v,r + v6t,t + 6t) = {Y,S{r-r,{t))6{v-v,{t)~ai{t))) 

i 

= - a, ■ d + ^aia, : dd + . . .]5{x - Xi{t))) (3) 

i 

where d = d/dv is a. gradient in the variable and (:) denotes a double contraction. As 
the total force is small for small 6t, the argument of the (5-function may be expanded in 
powers of aj. A subsequent average over the random noise gives, 

fiv, r + v5t, t + 5t)- fiv, r, t) = 5tJ,{f) + {5tfj^{f) + Omf) + ■ ■ • , (4) 

where the collision ieims Ji{f) and J^2{f) can be calculated straightforwardly. Here we only 
quote the dominant term of 0{6t) explicitly: 

J^{f) = d- J dx'k{x,x')f^^\x,x') + ]^dd : j dx'(T{x,x')(T{x,x')f'^^\x,x'). (5) 
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The term (/) contains 0{d^) with n = 2,3,4 and involves the higher order distribution 
functions Z^^-* and f^^\ Equation constitutes the first equation of the BBGKY-hierarchy 
for the DPD fluid with discrete time steps, relating the change in / to higher order distri- 
bution functions. 

To obtain a closed kinetic equation for f{x, t), the so-called the Fokker-Planck-Boltzmann 
(FPB) equation, we assume molecular chaos, i.e. f^'^\x,x',t) = f{x,t)f{x',t), etc. The 
dominant collision term contains the systematic frictional force as well as the random 
force. Next consider the propagation term, which is the Ihs of (^. In the limit of small St it 
simplifies to the usual streaming term dtf + v ■ V/, and the FPB equation becomes 

dtf{x) + vVf{x)=I{f). (6) 

Here V = 9/9r is a spatial gradient and the collision term becomes: 

/(/) = d- J dx'-f{x, x')f{x')f{x) + ^dd : J dx'(T{x, x')a{x, x')f{x')f{x). (7) 

We will concentrate on the case of continuous time. The most fundamental properties of 
the FPB equation are: it conserves particle number N = J dxf{x,t) and total momentum 
P = J dx vf{x,t), and obeys an if-theorem, dtJ^ < 0, for the total free energy: J-'{f) = 
J dx{\mv'^ + 6*0 In provided the detailed balance conditions 0, 6q = mcr^/27 and 
wd{t) = w1j{r) = w{r), are satisfied. The constant 6q depends only on the model parameters. 
In absence of conservative forces, the if-theorem implies the existence of a unique equilibrium 
solution, /o(a;) = UQipo^v) where ^po{v) = (m/27r6'o)'^''^ exp[— |mt>^/6'o] is a Maxwellian, and 
no = N/V is the density. If the detailed balance conditions are violated, the proof of the 
if-theorem breaks down. The modified algorithm ||^ satisfies the detailed balance conditions 
in the limit cit — > 0, but the original one []T[ does not. 

Our goal here is to derive the macroscopic evolution equations, describing the fluid dy- 
namics on large spatial and temporal scales, i.e. the Navier-Stokes equation, and to obtain 
explicit expressions for the thermodynamic and transport properties in terms of density n, 
temperature Oq, friction 7 and range Rq. 

As the interactions in the DPD fluid do not conserve energy, only mass density, p{r, t) = 
mn(r,t) = J dvf{r,v,t) , and momentum density p{r,t)u{r,t) = J dvf{r,v,t)v, satisfy 
local conservation laws, where u{r, t) is the local flow velocity of the fluid. The approach 
to equilibrium proceeds in two stages: first, a rapid relaxation {kinetic stage) to local equi- 
librium within a characteristic kinetic relaxation time to- The energy density e{r,t) decays 
within the same short kinetic stage to i^dO^nir, t), where ^0 is the global equilibrium temper- 
ature 9q = ksTo. In the subsequent hydrodynamic stage the time evolution of f{x,t) in the 
DPD fluid is only determined by the slow fields n{r,t) and u{r,t), i.e. f{v\n{r,t),u{r,t)). 
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but the temperature is thermostated at the global equilibrium value 60. All processes proceed 
isothermally. As a DPD fluid is not able to sustain a temperature gradient on hydrodynamic 
time scales, there is no heat current proportional to a temperature gradient and there is no 
heat conductivity. 

Using the Chapman- Enskog method / can be determined perturbatively, / = /o + 
/i/i + . . ., as an expansion in powers of a small parameter, /i ~ ^oV, which measures 
the variation of the macroscopic parameters over the characteristic kinetic length scale, 
io = vto oc v/'-yn with v = ^9o/m the mean velocity. The first term /o is the local equilib- 
rium distribution, /o = n{r,t) (m/(27r^o))'^^^ exp [— m(w — tt(r, t))^/(20o)] • The next term 
/i, linear in the gradients, is obtained as the solution of dtfo + v ■ V/o = (dlo/df) j^fi, 
where the time derivatives dtu and dtU are eliminated using the lowest order hydrodynamic 
equations (Euler equations). So, we first need to consider the local conservation laws. 

Integration of over v yields at once the continuity equation, dtn = —'V-nu. Similarly, 
we derive the conservation law for the momentum density by multiplying (^) by v and by 
integrating over v. The result after a partial v- integration is 

dtpu = —V- J dvvvf{x)—m J dvdx'j{x,x')f{x')f{x) 

= -V- (pHH + nK + no). (8) 

The kinetic part of the pressure tensor is the momentum flux in the local rest frame of the 
fluid IIk = / dvmvvf — puu = J dvmW f, where V = v — u{r, t) is the peculiar velocity. 
Substitution of / = /o + /i/i gives to lowest order IIk — "'^^'qI + Oi^p), where nO^ is the local 
equilibrium pressure. The next term on the rhs of (§) defines the dissipative part Hd of the 
momentum flux. It reduces to 

V ■ Hd = m7 J dRw{R)KR ■ (w(r) - u{r'))n{r)n{r'), (9) 

where r' = r — R. It is a typical collisional transfer contribution, resulting from the non- 
locality of the collision operator |p. Expansion of u{r') and n{r') around r shows that the 
rhs of (P) is 0{p^). It contributes to the viscosities but not to the Euler equations. After 
some algebra we obtain FId = —2riYyD — C^^SJ ■ nl, where the dissipative parts of the shear 
and bulk viscosity are identified as: 

r/D =mncJo(i?^)^/2(rf + 2) ; = mnujQ{R^)^/2d. (10) 

where we define the kinetic relaxation rate ujq = I/Iq = {n'-)/d)[w] with [w] = J dIlw{R) the 
effective volume of the action sphere, and where < >w= [R'^w]/[w^] is of order R^. 
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Following the standard Chapman- Enskog scheme, the solution /i is found to be propor- 
tional to Vw. The corresponding part of the kinetic pressure tensor then has the form 



Hk,! = J dvmVVfi = -2r]KD - CkV • ul, (11) 



where the rate of shear tensor D^^ is defined as the traceless symmetric part of VaUp. For 
the kinetic part of the viscosities we obtain the explicit results 

rjK = n9o/2uJo ; CK = n9o/duJo, (12) 

which are new results. The final transport coefficients are then r/ = r^o +^?k and C = Co + Ck- 
For further details we refer to 

The contribution is identical to the estimate for the total shear viscosity of the DPD 
fluid, calculated in |^ on the basis of the continuum approximation to the equations of motion 
for the DPD particles. Hoogerbrugge and Koelman have also shown that the viscosity found 
in numerical simulations does indeed approach riD for large values of wy. Here we have 
confirmed and rederived the formula for the viscosity of |]l| within a kinetic theory context. 

We conclude this paragraph by listing the complete results for the viscosities of the DPD 
fluid in the continuous case {6t —* 0), 



^ = >o 7T^ + - ; C = >or^ + - . (13) 





The traversal time of an action sphere is defined through = m(i?^)^/6'o = {R^)w/v'^. 
The expressions for the transport coefficients ([l^) involve the two intrinsic time scales of the 
DPD fluid: the characteristic kinetic time to = V'^O; and the traversal time which is of 
order Rq/v. In the parameter range > to, the dissipative viscosities ?7d and Co dominate 
and the estimate of |l|] is a reasonable one. In the range < to the kinetic viscosities r^K and 
Ck dominate. In a similar way the coefficient of self-diffusion D is obtained as -D = 9o/uJom. 

Results of numerical simulations for the kinematic viscosity C,/{nm) show |^ that the the- 
oretical results correctly describe the large and small ujq ~ n7 dependence, but a background 
contribution, roughly independent of ojq seems to be lacking in the theoretical predictions. 
There may be several reasons for this discrepancy, which are currently under investigation: 

(i) Effects of the discreteness of 5t, which alter the equilibrium temperature 0. 

(ii) Poor convergence of the Chapman- Enskog expansion, resulting from a poor separation 
of the "fast" kinetic and "slow" hydrodynamic timescales. 

(iii) The system size L, measured in units of the interaction range Rq, as used in the simula- 
tions of 10, may be too small and one is observing finite size effects related to wavenumber 
dependent transport coefficients ( generalized hydrodynamics). 
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(iv) Breakdown of the stosszahlansatz, caused by the small net momentum transfer in DPD 
colhsions. Consequently, dynamic correlations, built up by sequences of correlated binary 
("ring") collisions, may not be negligible. 

The main results of the present paper are the explicit predictions for the viscosities and 
self-diffusion coefficient of the DPD fluid in terms of the model parameters: density n, friction 
7, noise strength a or equivalently temperature 60 = mcr^/27 and range Rq. The expressions 
for rj and ( essentially depend on the two intrinsic time scales of DPD: the characteristic 
kinetic relaxation time to and the traversal time of the action sphere. 

The method used in eqs. (0-^ for deriving equations of motion for the reduced distribtion 
functions, can be directly applied to the A^-particle distribution function, P{yi, y2, ■ ■ ■ , Un, t) = 
(n.^^^5{yi — Xi{t))), and yields the discrete time analog of the A^-particle Fokker-Planck equa- 
tion for continuous time 0. Our equation for finite step size reduces to the Fokker-Planck 
equation of in the limit 6t 0. 
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